#include <iostream>
#include <limits>

double EPSILON = std::numeric_limits<double>::epsilon();

class FunctionBase
{
public:
    virtual double operator()(double x) = 0;

    virtual double differentiate(double x){
        double result;
        result = (this->operator()(x) - this->operator()(x - EPSILON)) / EPSILON;
        return result;
    }
};

template <class FuncType>
class NewtonInterpolator
{
private:
    double coefficients[100];    //保存差分系数
    double dataPoints[100];      //保存给定的数据点
    int pointCount;              //数据点的数量
    FuncType function;           //函数对象

public:
    NewtonInterpolator(double *points, int size, FuncType inputFunction) {
        function = inputFunction;
        pointCount = size;
        double tempCoeff[100][100];
        
        for(int index=0; index < size; index++){
            tempCoeff[index][0] = function(points[index]);
            dataPoints[index] = points[index];
        }
        
        for(int i = 1; i < size; i++){
            for(int j = 1; j <= i; j++){
                tempCoeff[i][j] = (tempCoeff[i][j-1] - tempCoeff[i-1][j-1]) / (points[i] - points[i-j]);
            }
        }
        
        for (int idx = 0; idx < size; idx++){
            coefficients[idx] = tempCoeff[idx][idx];
        }
    }


    double interpolate(double x) {//计算的差分系数来评估在点x处的牛顿插值多项式的值。
        double result = coefficients[pointCount-1];
        
        for(int i = 0; i < pointCount - 1; i++){
            result = result * (x - dataPoints[pointCount-i-2]) + coefficients[pointCount-i-2];
        }
        
        return result;
    }

    double *getCoefficients(){//返回差分系数数组
        return coefficients;
    }
};

class Newton_1
{
private:
    double coefficients[100];
    double dataPoints[100];
    int pointCount;
    double functionData[100];

public:
    Newton_1(double *points, int size, double *data) {
        pointCount = size;
        double tempCoeff[100][100];
        
        tempCoeff[0][0] = data[0];
        dataPoints[0] = points[0];

        for(int index=1; index<size; index++) {
            if(points[index] == points[index-1]) {
                tempCoeff[index][0] = tempCoeff[index-1][0];
            } else {
                tempCoeff[index][0] = data[index];
            }
            dataPoints[index] = points[index];
            functionData[index] = data[index];
        }

        int flag = 0;
        for(int i=1; i<size; i++) {
            flag = 0;
            for(int j=1; j<i+1; j++) {
                if(points[i] == points[i-j]) {
                    if(flag == 1) {
                        tempCoeff[i][j] = tempCoeff[i-1][j];
                    } else {
                        tempCoeff[i][j] = functionData[i];
                        flag = 1;
                    }
                } else {
                    tempCoeff[i][j] = (tempCoeff[i][j-1] - tempCoeff[i-1][j-1]) / (points[i] - points[i-j]);
                    flag = 0;
                }
            }
        }

        for(int idx=0; idx<size; idx++) {
            coefficients[idx] = tempCoeff[idx][idx];
        }
    }

    double interpolate(double x) {
        double result = coefficients[pointCount-1];
        for(int i = 0; i<pointCount-1; i++) {
            result = result * (x - dataPoints[pointCount-i-2]) + coefficients[pointCount-i-2];
        }
        return result;
    }

    double *getCoefficients() {
        return coefficients;
    }
};
